Climate change‐induced distributional change of medicinal and aromatic plants in the Nepal Himalaya

Abstract Medicinal and aromatic plants (MAPs) contribute to human well‐being via health and economic benefits. Nepal has recorded 2331 species of MAPs, of which around 300 species are currently under trade. Wild harvested MAPs in Nepal are under increasing pressure from overexploitation for trade and the effects of climate change and development. Despite some localized studies to examine the impact of climate change on MAPs, a consolidated understanding is lacking on how the distribution of major traded species of MAPs will change with future climate change. This study identifies the potential distribution of 29 species of MAPs in Nepal under current and future climate using an ensemble modeling and hotspot approach. Future climate change will reduce climatically suitable areas of two‐third of the studied species and decrease climatically suitable hotspots across elevation, physiography, ecoregions, federal states, and protected areas in Nepal. Reduction in climatically suitable areas for MAPs might have serious consequences for the livelihood of people that depend on the collection and trade of MAPs as well as Nepal's national economy. Therefore, it is imperative to consider the threats that future climate change may have on distribution of MAPs while designing protected areas and devising environmental conservation and climate adaptation policies.

surge of demand for natural health products and herbal drugs in recent times, the trading of MAPs is growing rapidly worldwide (Chen et al., 2016). In 2003, the annual global market for herbal medicines was estimated at US$60 billion, and by 2012, the global industry in Traditional Chinese Medicine (TCM) alone was reported to be worth US$83 billion (IPBES, 2019;Willis, 2017). MAPs are also a source of income for global rural populations through collection and sales after gathering from the uncultivated environments (Barata et al., 2016).
The Himalayan region, one of the world's biodiversity hotspots, has the highest concentration of MAP species (Kala, 2000;Olsen, 2005;Rai et al., 2000). One of the Himalayan countries, Nepal comprises 2331 recorded species of MAPs, of which around 300 species are currently under trade (Pyakurel et al., 2019;Rokaya et al., 2012). Likewise, two other Himalayan countries, China and India, have 11,146 and 7500 species of MAPs, respectively (Chen et al., 2016). The export of MAPs, including raw and processed plant products like Ayurvedic and traditional medicines produced from MAPs, in Nepal was worth around US$ 60.09 million in 2014, with an average annual export equivalent to 13.23 thousand tons (Ghimire et al., 2016). Although it shared <1% of global supplies, it contributed to approximately 5% of Nepal's GDP and 10% of the revenue collected from the forestry sector (Price, 2004). Along with a significant contribution to the national economy, MAPs also provide supplementary income and medicine for healthcare to rural households in Nepal (Larsen & Smith, 2004). In the mountainous regions of Nepal, commercial trade of wild alpine medicinal plants played an important role in the rural livelihoods, contributing on average 12% of annual household income (Olsen & Overgaard Larsen, 2003).
Harvesting of a highly valued medicinal species such as the caterpillar fungus (Ophiocordyceps sinensis) provides even much higher income to rural households; up to 65% of the total household income of the mountain communities in Nepal (Shrestha & Bawa, 2014).
In recent decades, medicinal plants are under increasing pressure from overexploitation for trade and the effects of climate change and development (Kling, 2016). This negatively affects the large portion of the global population who rely on natural medicines and reduces the potential to identify new medicinal compounds (Hopping et al., 2018;IPBES, 2019). In the Himalaya, climate change has impacted and will likely continue to impact biodiversity and ecosystems to various degrees (Bhattacharjee et al., 2017;Shrestha et al., 2012;Xu et al., 2009). Changes in community composition, distributional range, and growth pattern of a few species, including medicinal plants, were reported or predicted in Nepal due to climate change. For example, tree lines in the highaltitude region of the Himalaya are shifting upward Lamsal, Kumar, Shabani, & Atreya, 2017;Tiwari & Jha, 2018), including the growth of vegetation in sub-nival areas of Nepal's Himalayan region (Anderson et al., 2020). The decline and increase of suitable habitats of two medicinal plant species, namely Fritillaria cirrhosa and Lilium nepalense, have been predicted (Rana et al., 2017). Conversely, the suitable habitat of the medicinal fungus, Ophiocordyceps sinensis, has been expected to expand in the future with climate change in Nepal (Shrestha & Bawa, 2014). A recent study showed that potentially suitable habitats of Dactylorhiza hatagirea, Paris polyphylla, and Taxus spp. will expand particularly toward the north of Nepal under future climate (Kunwar et al., 2021), whereas Rana et al. (2020)  This study enhances the knowledge and understanding of the distribution of 29 species of MAPs in Nepal under current and future climate using an ensemble of species distribution models. We identified the current climatic envelope and estimated the future distribution of 29 species of MAPs. We further examined how the change in the distribution of medicinal plant hotspots (areas with a suitable climatic niche for the maximum number of species superimposed) will occur according to elevation, physiography, ecoregions, federal states, and protected areas. We also discussed how the future distributional change of the investigated species would affect the supply of medicinal raw materials to local people and the export industry.
The results of this study will be crucial to devise conservation strategies for MAPs in Nepal, especially at a time when the conservation of MAPs from overexploitation and climate change is pertinent.

| Study area
The study area covers the entire country of Nepal that lies at the center of the Himalaya biodiversity hotspot covering the area of 147,181 km 2 ( Figure 1). The country is divided into five physiographic zones, has sub-tropical to alpine climates and elevation ranges from 64 to 8848 m-the Mount Everest, seven federal states, and nine Global 200 ecoregions (Olson et al., 2001). In addition, 24% of the country's land area is covered by protected areas that comprise 12 national parks, one wildlife reserve, one hunting reserve, six conservation areas, and 13 buffer zones. The country is rich in biodiversity harboring little more than 6000 species of flowering plants of which 312 are endemic to Nepal (Tiwari et al., 2019).

| Species description and occurrence records
Twenty-nine medicinal and aromatic plant species found in Nepal Himalaya were selected based on their wide medicinal usage, conservation status, trading value, and availability of occurrence data ( Table 1) for Zanthoxylum armatum). Utilizing Google Earth, we geocoded 30 occurrence records that only have the names of collection localities. Data collected from the various sources were compiled, and duplicates and dubious records (e.g., records that fall 1000 m out from the reported elevation range of the species) were removed. Survey biases often displayed by species distributional data could have implications for predicting species occurrence under changing environmental conditions (Dormann, 2007). Spatial autocorrelation of sampling effort between training and test data inflates the prediction accuracy (Veloz, 2009). Therefore, spatial filtering is conducted to reduce sampling biases and model over-fitting (Boria et al., 2014;Dimson et al., 2019;Kramer-Schadt et al., 2013). Therefore, multiple presence locations in the same grid of ~1 km 2 spatial resolution (unit of analysis of this study) were removed and retained only one record per grid using the spatial filtering tool of SDMTOOLBOX 2.3 (Brown, 2014). Remaining 922 occurrence records were used in our ensemble modeling after removing erroneous and duplicated records.

| Environmental variables and model used
We downloaded 19 bioclimatic variables from the WorldClim data set (www.world clim.org) at 30 arc sec (~1 km 2 ) resolution. These bioclimatic variables were derived from monthly values of minimum, average and maximum temperature, and precipitation from 1970 to 2000 (Hijmans et al., 2005). The use of such relatively fine resolution of climate data is appropriate for regions with complicated topography, such as the Himalaya, where climatic conditions change significantly over a short distance. We analyzed a multicollinearity test among 19 bioclimatic variables and removed highly correlated variables (r > .70). Strong collinearity between the variables in predictive modeling could influence the overall model outcome by placing high emphasis on two or more highly correlated variables (Baldwin, 2009), resulting in misinterpretation.

| Ensemble modeling
An ensemble modeling of species distributions involves simulations across more than one set of initial conditions, model classes, model parameters, and boundary conditions (Araújo & New, 2007). The ensemble model accounts for uncertainties in predictions of different algorithms and uses a wide range of approaches to test models (Aguirre- Gutiérrez et al., 2017;Thuiller et al., 2009). However, a single-algorithm modeling method, MaxEnt can produce distribution maps of comparable accuracy to ensemble models (Kaky et al., 2020).
We used ensemble modeling because this consensus approach can often perform better than a single algorithm (Araújo & New, 2007;Thuiller et al., 2009). The analysis was conducted in R environment  (Hao et al., 2019).
Due to the unavailability of real absence data, we followed Barbet-Massin et al. (2012) and used 5000 pseudo-absences selected randomly for each repetition outside a buffer of 10 km from the presence points. The models were calibrated by using 70% of the occurrence points (presence and pseudo-absence) as training data and evaluated by using the remaining 30% as testing data (Araújo et al., 2005). We repeated the process of pseudo-absence generation three times by three evaluation runs per species, resulting in a total of 63 models per species (seven models, three evaluation runs and three pseudo-absence selection procedures) under each climate scenario. However, the use of pseudo-absence data might create inaccurate model performance (Liao & Chen, 2022). Therefore, the absence of real absence data is one of the limitations of this study.
We used True Skills Statistics (TSS) as an evaluation measure of model validation and predictive performance. TSS value ranges from −1 to +1 where +1 indicates a perfect agreement, and a TSS value below 0.4 indicates poor model discrimination (Allouche et al., 2006;Beaumont et al., 2016). Models with good predictive accuracy (TSS > 0.6) were used to build an ensemble from the projection outputs (Bellard et al., 2013;Gallien et al., 2012;Thuiller et al., 2009).
From the 63 individual models per species, we built ensemble models using a weighted-mean approach in which weights are awarded for each model proportionally to their evaluation metrics scores; hence, the discrimination is fair in this approach (Marmion et al., 2009).
Binary maps (suitable and unsuitable) were created using the optimal threshold that maximizes the TSS score as a cutoff value, which then converted the projected occurrence probabilities during the crossvalidation procedure (Allouche et al., 2006;Liu et al., 2013;Marmion et al., 2009). This threshold is unaffected by the prevalence of spe- net/c/world -datab ase-on-prote cted-areas) were used.

| RE SULTS
We evaluated the performance of models by TSS performance matrix ( Figure S1). The average TSS value of our models is 0.67 indicating good predictive accuracy. We excluded the models with TSS value <0.6 for building the ensemble-models.

| Suitable areas for individual species
On average, 5821 km 2 of area was predicted to be climatically suitable for individual species modeled. However, the results revealed a wide species-specific variation in climatically suitable areas ( Figure 2a). Zanthoxylum armatum was predicted to have the largest suitable area that is 45 times higher than the area suitable for Swertia chirayita, the species with the smallest suitable area ( Table 2). The other species with more than 10,000 km 2 of suitable area included

Valeriana hardwickei, Paris polyphylla, and Rhododendron anthopogon.
For eight species, including Swertia chirayita, Dactylorhiza hatagirea, Astilbe rivularis, and Rubia manjith, the suitable area was predicted to be <2000 km 2 . In the future, the average suitable area would decline by 10.4% with 5215 km 2 being predicted to be suitable for individual species. Compared to the suitable area predicted under current climate, there will be a decline in suitable areas of 19 species (65.5%) in the future with >50% reduction predicted for the spe- suitable area under future climate was predicted to be only marginal (−1.2%) for Nardostachys jatamansi (critically endangered and CITES II listed species) but significantly large for the remaining two species; 74% decline of Podophyllum hexandrum and 68% increase of Dioscorea deltoidea.

| Extent and hotspots
The extent of climatically suitable areas, i.e., suitable for at least one of the 29 MAP species studied, was predicted to be 57,306 km 2 under current climate with a 4% reduction in the future (Figure 2b).
About 29.5% of the current extent of suitable areas was predicted to be hotspots for the studied species that would decline to 18.5% in the future (Figure 2c). Considering the hotspot, about 40% of the current hotspot areas will be lost in the future due to climate change.
The results showed that there would be an overall decline in the extent of suitable areas as well as hotspots, but the decline in hotspot area will be larger than the decline in the extent of the suitable area.

| Hotspots in protected areas
About 52% of the current total hotspot area for the entire study area is located within 13 protected areas (Table 3). Among the protected areas, Annapurna Conservation Area was predicted to have the largest hotspot areas (2813 km 2 ), followed by Gaurishanker Conservation Area (1360 km 2 ) and Langtang National Park (1008 km 2 ). In the future, the percentage of hotspot area within these protected areas would slightly increase and reach to about 60%. However, the hotspot area would decline in the majority (10 out of 13) of the protected areas.
Among the protected areas with the large hotspot areas, the highest decline was predicted in Annapurna Conservation Area (58%).

| Hotspots in global ecoregions
Across the global ecoregions, the hotspots of suitable areas are mainly concentrated in four ecoregions, namely, Eastern Himalayan alpine shrub and meadows, Eastern Himalayan subalpine conifer forests, Western Himalayan alpine shrub and meadows, and Eastern Himalayan broadleaf forests ( Table 4). Each of these ecoregions had areas >2000 km 2 that were predicted to be hotspots. In the future, hotspot area would decline in majority (7 out of 9) of the ecoregions with the highest decline occur in the Western Himalayan alpine shrub and meadows.

| Hotspots along elevation gradients
Distribution of hotspot areas varied significantly with elevation ( Figure 3). The hotspot areas were mostly concentrated between 2000 and 5000 m a.s.l. with the largest hotspot area being predicted between 4000 and 4500 m a.s.l. For the studied species, the hotspot was not predicted below 1000 m a.s.l. In the future, a decline in hotspot area was predicted for the elevation range between 2000 and 5500 m a.s.l. and an increase in hotspot area below and above that elevation range.
The hotspot areas are mainly concentrated in High Mountain regions, followed by the Middle Mountains and Hill. The High Mountain region would have the highest proportion of hotspot area but would have a 47% reduction of hotspot area under future climate scenario (Figure 4a).

| Hotspots across federal states
There was also significant variation in the distribution of hotspot areas across the federal states ( Figure 4b). The highest percentage of hotspots was predicted to lie within Gandaki province (30.6%), followed by Bagmati province (27.5%) and Province 1 (23.8%). For the studied species, there was no hotspot area in Madhesh province (Province 2) and very little hotspot area in Lumbini province (Province 5). All seven provinces would lose hotspot areas with the highest loss in Lumbini province (97.6%) and the lowest loss in Bagmati province (12.2%).

| D ISCUSS I ON S AND CON CLUS I ON S
To our knowledge, we presented a comprehensive analysis of climat- Nepal has already experienced warming, increased annual precipitation, and increased extreme events, including an increased number of hot days and nights and heavy precipitation . By the end of 21 st century, Nepal could face an increased mean annual temperature between 1.7-3.6°C and 11%-23% increased precipitation (MoFE, 2019). Climate change in Nepal has already impacted landscape phenology , the distribution of forests (Thapa et al., 2016) and tree line positions (Tiwari & Jha, 2018 (Ghimire et al., 2008). These excessive and detrimental harvesting techniques negatively impact the population of several species of medicinal plants (Ghimire et al., 2008;Shrestha & Bawa, 2013).
Additionally, extinction risk for herbaceous plants is higher than for woody plants (Yan & Tang, 2019). However, the extinction risk depends on several factors such as harvesting pressure and parts harvested. Therefore, the existing detrimental harvesting practices combined with the predicted decline of climatically suitable areas of MAPs will pose additional threats to the MAPs population and associated trade in Nepal. This decline could affect exports of MAP species as MAPs are currently exported to 50 countries from Nepal (Ghimire et al., 2016). Furthermore, Gairola et al. (2010) and Gupta and Chaturvedi (2019)  Protected areas should be able to maintain a long-term dynamics of biodiversity change (Pressey et al., 2007). The climatically suitable hotspots of MAPs in nine protected areas will decline under future climate change scenarios, indicating that the coverage of the existing protected areas might not be suitable to conserve highly traded MAP species in the future. Furthermore, the existing protected areas are not fully representative and failed to incorporate diverse topography, ecosystems, vegetation, flora, and fauna of the country (Shrestha et al., 2010). Therefore, ensuring adequate representation of topography, ecosystems and species is essential to improve the existing protected areas system of Nepal. Our study facilitates future conservation planning of the existing conservation areas by identifying shifts in MAPs distribution within and outside protected areas under climate change.
Our results showed that future climate change will reduce climatically suitable areas for majority of the traded MAPs in Nepal.
The MAPs in Nepal have been a major contributor to traditional health care, household income, and export. Excessive and destructive harvesting practices have raised a concern toward the conservation of various species of MAPs particularly herbaceous perennial species that already have higher risk of extinction. The Government of Nepal has already listed some MAP species under the national conservation/protection list that imposes a ban on collection, trade, and export. Our results showed that climatically suitable areas of the majority of traded MAP species will be reduced with future climate change. Reduction in climatically suit-

ACK N OWLED G M ENTS
The data used in this study was collected by SKG in his various fieldworks that are supported by People and Plants Initiative, WWF Nepal, WWF UK, ICIMOD, GIZ, National Geographic Society, and Missouri Botanical Garden. UBS also thanks to National Geographic Society (grant number NGS-62058R-19) to support his work. We thank Nirmala Phuyal for providing occurrences of Zanthoxylum armatum and National Herbarium and Plant Laboratories Godavari (KATH) and Tribhuvan University Central Herbarium (TUCH) for allowing us to collect occurrence data. We are grateful to Lily Clark for English proofreading.

DATA AVA I L A B I L I T Y S TAT E M E N T
Upon acceptance, the data that support the findings of this study are openly available in Dryad.